Purpose: To measure the growth rate of all strains, starting with cultures grown to mid log phase
Protocol:
library(tidyverse)
library(xlsx)
library(broom)
library(rstatix)
library(ape)
library(formattable)
library(corrplot)
library(ggbeeswarm)
library(cowplot)
library(ggpubr)
library(grid)
library(gridGraphics)
`%!in%` <- negate(`%in%`)
Short cuts and file paths:
strains <- c("14018", "315-A", "C0011E4", "JCP7275", "C0084H9", "JCP8108", "JCP8017B", "C0040C2" , "UM35", "JCP8066", "C0093B3", "AMD", "C0096A1", "C0102A2", "UM224", "Gv5-1", "C0179B4", "C0056B5", "Gv101", "C0100B2", "C0101A1", "CMW7778B")
strains0 <- c("14018", "315-A", "C0011E4", "JCP7275", "C0084H9", "JCP8108", "JCP8017B", "C0040C2" , "UM35", "JCP8066", "C0093B3", "AMD", "C0096A1", "C0102A2", "UM224", "5-1", "C0179B4", "C0056B5", "101", "C0100B2", "C0101A1", "CMW7778B")
species <- c("Gardnerella vaginalis", "Gardnerella sp. 2", "Gardnerella sp. 3", "Gardnerella piotii", "Gardnerella leopoldii", "Gardnerella swidsinskii", "Gardnerella sp. 7", "Gardnerella sp. 8", "Gardnerella sp. 10", "Gardnerella sp. 11", "Gardnerella sp. 12")
batches <- c("1A1", "1A2", "1B1", "1B2", "2A1", "2A2", "3A1", "3A2")
figureOut <- "../../experiments_figures"
suppTableOut <- "../../experiments_figures/supplementary_posthoc_tables.xlsx"
Import data
midlog_raw_data <- read_csv("./data/gc_midlog_raw_data.csv")
curve_raw_data <- read_csv("./data/gc_curve_raw_data.csv")
Plot settings
theme_set(theme_bw())
speciesColor <- scale_color_manual(values=c("#CC79A7", "#D55E00", "#56B4E9", "#0072B2", "#E69F00", "#009E73",
"darkgreen", "gold3", "darkred", "darkblue",
"mediumpurple4"))
speciesFill <- scale_fill_manual(values=c("#CC79A7", "#D55E00", "#56B4E9", "#0072B2", "#E69F00", "#009E73",
"darkgreen", "gold3", "darkred", "darkblue",
"mediumpurple4"))
Propagate twice and then incubate to mid log. 1:10 Dilution for mid Log: Protocol 1: 4 hour incubation Protocol 2: 7 hour incubation Protocol 3: 15 hour incubation ## Clean data
midlog_df <- midlog_raw_data %>%
group_by(Number, BioRep, Species, Strain, Batch, Step) %>%
mutate(OD=OD-mean(.$OD[.$Number=="Blank"]),
OD=case_when(OD>0~OD,
OD<=0~0)) %>%
filter(Number!="Blank") %>%
mutate(Protocol=str_sub(Batch, 1, 1),
BatchNum=str_sub(Batch, 3,3),
Step=factor(Step, levels=c("Overnight 1", "Overnight 2", "Dilution Start", "Mid Log")),
Strain=factor(Strain, levels = strains),
Species=factor(Species, levels = species)) %>%
with_groups(c("Number", "BioRep", "Species", "Strain", "Batch", "Step", "Protocol", "BatchNum"), summarize, ODm=mean(OD), ODsd=sd(OD)) %>%
filter(Step!="Dilution Start") %>%
mutate(StrainN=as.numeric(Strain),
rectFill=case_when(StrainN %% 2 == 0 ~ "A",
StrainN %% 2 != 0 ~ "B")) %>%
ungroup %>%
mutate(Sample=paste(Strain, Batch, BioRep, sep="."))
Assess mid log growth for removing poor growers with OD600 < 0.05
midlog_df %>%
filter(Step=="Mid Log") %>%
ggplot(aes(x=Strain, y=ODm, shape=BioRep, color=Batch)) +
geom_point() +
geom_hline(aes(yintercept=0.05), linetype=2, color="gray") +
theme(axis.text.x = element_text(angle=45, hjust=1)) +
labs(x=NULL, y=bquote(OD[600]), shape="Biological Replicate", color="Experiment Batch")
Create quality filter for OD600 > 0.05
midLogFail <- midlog_df %>%
filter(Step=="Mid Log",
ODm < 0.05) %>%
.$Sample
midLogFail
## [1] "C0100B2.1A1.A" "C0100B2.1A2.A" "C0100B2.1A1.B" "C0100B2.1A2.B"
## [5] "C0096A1.2A2.B" "JCP8066.1A2.A" "JCP8066.1A2.B" "UM224.1A1.B"
## [9] "C0100B2.1B1.A" "C0100B2.1B1.B"
Apply quality filter and plot set-up growth endpoints
# main plot with quality filtering
overnightPlot <- midlog_df %>%
filter(Sample %!in% midLogFail) %>%
ggplot(aes(x=Strain, y=ODm)) +
geom_rect(aes(xmin=StrainN-.5, xmax=StrainN+.5, ymin = 0, ymax = 1, fill=rectFill), alpha = 0.2, show.legend = FALSE) +
geom_pointrange(aes(ymin=(ODm-ODsd), ymax=(ODm+ODsd), color=Species), position = position_quasirandom(), size=0.3, alpha=0.7, show.legend = FALSE) +
speciesColor +
scale_fill_manual(values=c("white", "gray"), labels=c("A","B")) +
ylim(0,1) +
facet_grid(.~Step) +
coord_flip() +
scale_x_discrete(limits=rev) +
theme(axis.text.x = element_text(size=9),
axis.text.y = element_text(size=8.5),
axis.title = element_text(size=13),
strip.text = element_text(size=13),
legend.text = element_text(size=11),
panel.grid = element_blank(),
legend.position = "none") +
labs(x=NULL, y=bquote(OD[600]), shape="Experiment Batch")
#annotation bars
protocolBar <- midlog_df %>%
select(Strain, Species, Protocol) %>%
unique %>%
ggplot() +
geom_bar(aes(x=Strain, y=2, fill=Protocol), stat="identity", width = 1.1) +
scale_x_discrete(limits=rev) +
scale_fill_manual(values=c("cornflowerblue", "burlywood1", "forestgreen")) +
coord_flip() +
theme_void() +
theme(legend.position = "none")
speciesBar <- midlog_df %>%
select(Strain, Species, Protocol) %>%
unique %>%
ggplot() +
geom_bar(aes(x=Strain, y=1, fill=Species), stat="identity", width=1.1) +
scale_x_discrete(limits=rev) +
speciesFill +
coord_flip() +
theme_void() +
theme(legend.position = "none")
protocol_legend <- ggpubr::get_legend(protocolBar +
guides(color=guide_legend(ncol = 1)) +
theme(legend.position = "bottom",
legend.direction = "horizontal"))
species_legend <- ggpubr::get_legend(speciesBar +
guides(color=guide_legend(ncol = 1)) +
theme(legend.position = "bottom",
legend.direction = "horizontal"))
a <- plot_grid(overnightPlot, speciesBar, protocolBar, nrow = 1, align = "h", axis = "bt", rel_widths = c(1, 0.02, 0.02))
b <- ggdraw(plot_grid(a, as_ggplot(species_legend), nrow=2, rel_heights = c(1, 0.3)))
setUpGrowthFigure <- cowplot::plot_grid(b, as_ggplot(protocol_legend), ncol=1, nrow=2, rel_heights = c(1, 0.1))
setUpGrowthFigure
#ggsave(file.path(figureOut, paste(Sys.Date(), "Figure2_curves_setup_growth.png", sep="_")))
blanks0 <- curve_raw_data %>%
filter(Number=="Blank") %>%
select(well, Number, Time, OD, Batch)
blanks0 %>%
ggplot(aes(x=Time, y=OD, color=well)) +
geom_point() +
facet_wrap(~Batch, ncol=2, scales = "free")
Remove wells B12 and C12 from experiment batch 1A1
blanks0 %>%
mutate(expwell=paste(Batch, well, sep=".")) %>%
filter(expwell!="1A1.B12", expwell!="1A1.C12") %>%
ggplot(aes(x=Time, y=OD, color=well)) +
geom_point() +
facet_wrap(~Batch, ncol=2, scales = "free_x")
Assess blank variability
blankSummaryTime <- blanks0 %>%
mutate(expwell=paste(Batch, well, sep=".")) %>%
filter(expwell!="1A1.B12",
expwell!="1A1.C12") %>%
with_groups(c(Time, Batch), summarize, timeRange=max(OD)-min(OD))
blankSummaryWell <- blanks0 %>%
mutate(expwell=paste(Batch, well, sep=".")) %>%
filter(expwell!="1A1.B12",
expwell!="1A1.C12") %>%
with_groups(expwell, summarize, wellRange=max(OD)-min(OD))
timeBlankPlot <- blankSummaryTime %>%
ggplot(aes(x="Range per Timepoint", y=timeRange)) +
geom_violin(alpha=0) +
geom_quasirandom(alpha=0.5) +
stat_summary(fun = "mean", geom="crossbar", width=0.8, color="red", linetype=2) +
scale_y_continuous(breaks = seq(0, 0.0125, 0.0025), limits = c(0, 0.0125)) +
labs(x=NULL, y="OD Range")
wellBlankPlot <- blankSummaryWell %>%
ggplot(aes(x="Range per Well", y=wellRange)) +
geom_violin(alpha=0) +
geom_quasirandom(alpha=0.5) +
scale_y_continuous(breaks = seq(0, 0.0125, 0.0025), limits = c(0, 0.0125)) +
stat_summary(fun = "mean", geom="crossbar", width=0.8, color="red", linetype=2) +
theme(axis.text.y = element_blank()) +
labs(x=NULL, y=NULL)
plot_grid(timeBlankPlot, wellBlankPlot)
summary(blankSummaryTime$timeRange)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000000 0.002000 0.003000 0.003039 0.004000 0.009000
Limit of detection appears to be about 0.003.
blanks <- curve_raw_data %>%
mutate(expwell=paste(Batch, well, sep=".")) %>%
filter(Number=="Blank",
expwell!="1A1.B12",
expwell!="1A1.C12") %>%
with_groups(c(Batch, Time), summarize, blankOD=mean(OD))
curve <- curve_raw_data %>%
filter(Number!="Blank") %>%
left_join(blanks, by=c("Time", "Batch")) %>%
mutate(OD=OD-blankOD) %>%
select(-blankOD) %>%
mutate(Sample=paste(Strain, Batch, BioRep, sep="."))
curve %>%
filter(Sample %!in% midLogFail) %>% #filtering based on mid log failure
ggplot(aes(x=Time, y=OD, color=Species, shape=BioRep, linetype=BioRep)) +
geom_point(alpha=0.5) +
scale_x_continuous(breaks = seq(0,48,12)) +
scale_y_log10() +
facet_wrap(~Strain) +
theme(axis.text = element_text(size=9),
axis.title = element_text(size=11),
strip.text = element_text(size=10),
legend.title = element_text(size=11)) +
labs(x="Time (Hours)", y=bquote(OD[600]), color="Species", shape = "Biological Replicate", linetype="Biological Replicate")
In addition to samples removed for quality control at mid log step:
1A1: Remove well E9 (one technical replicate for 14018 biological replicate B)
1B1: Remove JCP8108 biological replicate A
curveFilt <- curve %>%
mutate(ExpWell=paste(Batch, well, sep=".")) %>%
filter(ExpWell!="1A1.E9",
Sample %!in% c(midLogFail, "JCP8108.1B1.A"))
curveMeans <- curveFilt %>%
filter(OD>0) %>% # remove any technical replicates below detection
with_groups(c(Sample, Number, BioRep, Strain, Species, Time, Batch), summarize, mOD=mean(OD), sd=sd(OD)) %>%
with_groups(c(Sample, Batch), nest, Time=Time, mOD=mOD, sd=sd)
Use rolling regression to calculate growth curves
# create the rolling regression function
roll_regress <- function(x){
temp <- data.frame(x)
mod <- lm(temp)
temp <- data.frame(slope = coef(mod)[[2]],
slope_lwr = confint(mod)[2, ][[1]],
slope_upr = confint(mod)[2, ][[2]],
intercept = coef(mod)[[1]],
rsq = summary(mod)$r.squared, stringsAsFactors = FALSE)
return(temp)
}
rollRegressData <- curveFilt %>%
mutate(SampWell=paste(Sample, well, sep=".")) %>%
arrange(SampWell, Time) %>%
mutate(Strain=factor(Strain, levels=strains),
Species=factor(Species, levels=species))
First look at curves with smallest log phases
curveFilt %>%
filter(Strain %in% c("14018", "315-A", "C0011E4", "C0093B3")) %>%
ggplot(aes(x=Time, y=OD)) +
geom_point(alpha=0.5) +
facet_wrap(~Strain, scales = "free_x") +
scale_y_log10()
Choose 2 hours for rolling window size. Choose 0.006 as minumum filtering threshold.
# define window - here every 2 hours with measurements every 0.5 hours
num_points <- ceiling(2*60/(60*0.5))
num_points
## [1] 4
# run rolling regression on lnOD ~ time
models <- rollRegressData %>%
filter(OD>0.006) %>%
mutate(ln_od=log(OD)) %>%
group_by(SampWell) %>%
do(cbind(model = select(., ln_od, Time) %>%
zoo::rollapplyr(width = num_points, roll_regress, by.column = FALSE, fill = NA, align = 'center'),
Time = select(., Time),
ln_od = select(., ln_od))) %>%
rename_all(., gsub, pattern = 'model.', replacement = '')
# create predictions
preds <- models %>%
filter(., !is.na(slope)) %>%
group_by(SampWell, Time) %>%
do(data.frame(time2 = c(.$Time - 2, .$Time + 2))) %>%
left_join(., models) %>%
mutate(pred = (slope*time2) + intercept)
# get growth rates
growth_rate <- models %>%
group_by(SampWell) %>%
filter(slope == max(slope, na.rm = TRUE))
rollRegressPlotSubset <- rollRegressData %>%
group_by(Strain) %>%
sample_n(1, replace = FALSE) %>%
mutate(Strain=factor(Strain, levels=strains)) %>%
arrange(Strain) %>%
.$SampWell %>%
unique
# plot rolling regression
rollRegressData %>%
filter(OD>0.006,
SampWell %in% rollRegressPlotSubset) %>%
mutate(ln_od=log(OD),
Strain=factor(Strain, levels=strains),
SampWell=factor(SampWell, levels=rollRegressPlotSubset)) %>%
ggplot(aes(x=Time, y=ln_od)) +
geom_point(size=1) +
geom_line(aes(time2, pred, group = Time), col = 'red', (subset(preds, SampWell %in% rollRegressPlotSubset) %>% mutate(SampWell=factor(SampWell, levels=rollRegressPlotSubset)))) +
geom_segment(aes(x = Time, y = ln_od-3, xend = Time, yend = ln_od), (subset(growth_rate, SampWell %in% rollRegressPlotSubset) %>% mutate(SampWell=factor(SampWell, levels=rollRegressPlotSubset)))) +
geom_segment(aes(x = ln_od-3, y = ln_od, xend = Time, yend = ln_od), (subset(growth_rate, SampWell %in% rollRegressPlotSubset) %>% mutate(SampWell=factor(SampWell, levels=rollRegressPlotSubset)))) +
geom_text(aes(label= Strain, x=16, y=1.5), size=3, hjust="center", check_overlap = TRUE) +
facet_wrap(~SampWell) +
ylim(-8.5, 2.5) +
theme(strip.text = element_blank(),
strip.background = element_blank(),
axis.text = element_blank(),
axis.ticks = element_blank(),
panel.grid = element_blank(),
panel.border = element_rect(color = "black", size=1)) +
labs(x = NULL, y = NULL) -> growthCurveFigures
growthCurveFigures
ratePlot <- growth_rate %>%
separate(SampWell, c("Strain", "Batch", "BioRep", "Well"), sep="\\.", remove=FALSE) %>%
mutate(Sample=paste(Strain, Batch, BioRep, sep=".")) %>%
left_join(unique(curve[,c("Strain", "Species")])) %>%
with_groups(c(Sample, Species, Strain, Batch, BioRep), summarize, mSlope=mean(slope), sdSlope=sd(slope)) %>%
mutate(Strain=factor(Strain, levels=strains),
Species=factor(Species, levels=species)) %>%
ungroup %>%
ggplot(aes(x=Strain, y=mSlope)) +
geom_point(alpha=0) +
geom_boxplot(position = position_dodge2(width=1), alpha=0, color="black", show.legend = FALSE) +
geom_pointrange(aes(x=Strain, y=mSlope, ymin=mSlope-sdSlope, ymax=mSlope+sdSlope, color=Species), size=0.5, alpha=0.6, position = position_dodge2(width=1), show.legend = FALSE) +
ylim(0,1.5) +
speciesColor +
theme(axis.text.x = element_text(angle=45, hjust=1),
axis.title.y = element_text(size=10),
panel.grid = element_blank(),
legend.position = "none") +
labs(x=NULL, y="Growth Rate")
protocolBar0 <- midlog_df %>%
select(Strain, Species, Protocol) %>%
unique %>%
ggplot() +
geom_bar(aes(x=Strain, y=2, fill=Protocol), stat="identity", width = 1.1) +
scale_fill_manual(values=c("cornflowerblue", "burlywood1", "forestgreen")) +
theme_void() +
theme(legend.position = "none")
speciesBar0 <- midlog_df %>%
select(Strain, Species, Protocol) %>%
unique %>%
ggplot() +
geom_bar(aes(x=Strain, y=1, fill=Species), stat="identity", width=1.1) +
speciesFill +
theme_void() +
theme(legend.position = "none")
#Legends
protocol_legend <- cowplot::get_legend(
protocolBar0 +
guides(color = guide_legend(nrow = 1)) +
theme(legend.position = "bottom"))
species_legend <- cowplot::get_legend(
speciesBar0 +
guides(color=guide_legend(ncol = 1)) +
theme(legend.position = "bottom"))
#plots
plots <- plot_grid(protocolBar0, speciesBar0, ratePlot, ncol = 1, rel_heights = c(0.08, 0.08, 1), align = "v", axis = "bt")
plot_grid(protocol_legend, species_legend, plots, ncol=1, rel_heights = c(0.1, 0.25, 1))
qqplot <- growth_rate %>%
separate(SampWell, c("Strain", "Batch", "BioRep", "Well"), sep="\\.", remove=FALSE) %>%
mutate(Sample=paste(Strain, Batch, BioRep, sep=".")) %>%
left_join(unique(curve[,c("Strain", "Species")])) %>%
with_groups(c(Sample, Species, Strain, Batch, BioRep), summarize, mSlope=mean(slope), sdSlope=sd(slope)) %>%
ggplot(aes(sample=mSlope)) +
geom_qq_line() +
geom_qq()
qqplot
strain_growth_rate_DF <- growth_rate %>%
separate(SampWell, c("Strain", "Batch", "BioRep", "Well"), sep="\\.", remove=FALSE) %>%
mutate(Sample=paste(Strain, Batch, BioRep, sep=".")) %>%
left_join(unique(curve[,c("Strain", "Species")])) %>%
with_groups(c(Sample, Species, Strain, Batch, BioRep), summarize, mSlope=mean(slope), sdSlope=sd(slope)) %>%
mutate(Strain=factor(Strain, levels=strains))
strain_growth_rate_DF %>%
anova_test(mSlope~Strain) %>%
formattable()
| Effect | DFn | DFd | F | p | p<.05 | ges |
|---|---|---|---|---|---|---|
| Strain | 21 | 61 | 14.494 | 1.79e-16 |
|
0.833 |
growthRatePostHoc <- strain_growth_rate_DF %>%
tukey_hsd(mSlope~Strain, p.adjust.method = "BH")
# growthRatePostHoc %>%
# arrange(p.adj) %>%
# write.xlsx(file=suppTableOut, sheetName="TableS3_growth_rate_posthoc", col.names=TRUE, append=FALSE)
min(growthRatePostHoc$p.adj)
pvals <- growthRatePostHoc %>%
select(group1, group2, p.adj) %>%
dplyr::rename(var1=group2, var2=group1, cor=p.adj) %>%
cor_spread() %>%
column_to_rownames(var="rowname") %>%
as.matrix()
growthRatePostHoc %>%
select(group1, group2, estimate) %>%
dplyr::rename(var1=group2, var2=group1, cor=estimate) %>%
cor_spread() %>%
column_to_rownames(var="rowname") %>%
as.matrix() %>%
corrplot(is.corr=FALSE, type = "lower", tl.col="black", p.mat = pvals, insig = "label_sig", sig.level = c(0.001, 0.01, 0.05), pch.cex = 0.9)
#bg="light gray", cl.pos="n",
grid.echo()
growthRatePostHocPlot <- grid.grab()
Load tree
We choose the weights as \(_{wij} = 1/d_{ij}\) , where the d’s is the distances measured on the tree:
We can now perform the analysis with Moran’s I:
carryingCapacity <- curveFilt %>%
with_groups(c("well", "Number", "BioRep", "Species", "Strain", "Batch", "Sample", "ExpWell"), summarize, OD=max(OD)) %>%
with_groups(c("Number", "BioRep", "Species", "Strain", "Batch", "Sample"), summarize, maxOD=mean(OD), maxODsd=sd(OD)) %>%
mutate(Strain=factor(Strain, levels=strains),
Species=factor(Species, levels=species))
carryingCapacityPlot <- carryingCapacity %>%
mutate(Strain=factor(Strain, levels=strains)) %>%
ggplot(aes(x=Strain, y=maxOD)) +
geom_boxplot(alpha=0) +
geom_pointrange(aes(ymin=maxOD-maxODsd, ymax=maxOD+maxODsd, color=Species), size=0.5, alpha=0.6, position = position_dodge2(width=1), show.legend = FALSE) +
speciesColor +
theme(axis.text.x = element_text(angle=45, hjust=1),
axis.title.y = element_text(size=10),
panel.grid = element_blank(),
legend.position = "none") +
labs(x=NULL, y=bquote("Carrying Capacity (OD"[600]~")"))
plots2 <- plot_grid(protocolBar0, speciesBar0, carryingCapacityPlot, ncol = 1, rel_heights = c(0.08, 0.08, 1), align = "v", axis = "bt")
plot_grid(protocol_legend, species_legend, plots2, ncol=1, rel_heights = c(0.1, 0.25, 1))
qqplot <- carryingCapacity %>%
ggplot(aes(sample=maxOD)) +
geom_qq_line() +
geom_qq()
qqplot
strain_growth_rate_DF <- carryingCapacity %>%
mutate(Strain=factor(Strain, levels=strains))
strain_growth_rate_DF %>%
anova_test(maxOD~Strain) %>%
formattable()
| Effect | DFn | DFd | F | p | p<.05 | ges |
|---|---|---|---|---|---|---|
| Strain | 21 | 61 | 6.81 | 1.96e-09 |
|
0.701 |
carCapacPostHoc <- carryingCapacity %>%
tukey_hsd(maxOD~Strain, p.adjust.method = "BH")
# carCapacPostHoc %>%
# arrange(p.adj) %>%
# write.xlsx(file=suppTableOut, sheetName="TableS4_car_capac_posthoc", col.names=TRUE, append=TRUE)
min(carCapacPostHoc$p.adj)
pvals <- carCapacPostHoc %>%
select(group1, group2, p.adj) %>%
dplyr::rename(var1=group2, var2=group1, cor=p.adj) %>%
cor_spread() %>%
column_to_rownames(var="rowname") %>%
as.matrix()
carCapacPostHoc %>%
select(group1, group2, estimate) %>%
dplyr::rename(var1=group2, var2=group1, cor=estimate) %>%
cor_spread() %>%
column_to_rownames(var="rowname") %>%
as.matrix() %>%
corrplot(is.corr=FALSE, type = "lower", tl.col="black", p.mat = pvals, insig = "label_sig", sig.level = c(0.001, 0.01, 0.05), pch.cex = 0.9)
grid.echo()
carCapacPostHocPlot <- grid.grab()
We can now perform the analysis with Moran’s I:
figure3B <- plot_grid(protocolBar0, speciesBar0, ratePlot + theme(axis.text.x = element_blank(), axis.ticks.x = element_blank()), carryingCapacityPlot, ncol = 1, rel_heights = c(0.08, 0.08, .75, 1), align = "v", axis = "bt")
plot_grid(ggplot()+geom_blank()+theme_void(), growthCurveFigures, figure3B, rel_widths = c(0.08, 1.2, 1), nrow = 1, labels = c("A", "", "B"), label_size = 15)
#ggsave(file.path(figureOut, paste(Sys.Date(), "Figure3_growth_rates.png", sep = "_")))
plot_grid(ggplot()+geom_blank(), growthRatePostHocPlot, ggplot()+geom_blank(), carCapacPostHocPlot, nrow = 1, rel_widths = c(0.1, 1, 0.1, 1), labels = c(" ", "A", " ", "B"), label_size = 20)
#ggsave(file.path(figureOut, paste(Sys.Date(), "FigureS2_growth_posthocs.png", sep = "_")))
cfuTable1A1 <- tribble(
~Number, ~BioRep, ~Dilution, ~Colonies,
"1", "A", 1e-4, 34,
"1", "B", 1e-4, 132,
"2", "A", 1e-4, 70,
"2", "B", 1e-3, 166,
"3", "A", 1e-3, 83,
"3", "B", 1e-2, 242,
"4", "A", 1e-4, 32,
"4", "B", 1e-3, 31,
"5", "A", 1e-4, 48,
"5", "B", 1e-4, 23,
"6", "A", 1e-4, 76,
"6", "B", NA, NA,
"7", "A", 1e-3, 110,
"7", "B", 1e-4, 71,
"8", "A", 1e-3, 76,
"8", "B", 1e-5, 25)
cfuTable1A2 <- tribble(
~Number, ~BioRep, ~Dilution, ~Colonies,
"1", "A", 1e-3, 53,
"1", "B", 1e-4, 51,
"2", "A", 1e-3, 29,
"2", "B", 1e-3, 103,
"3", "A", 1e1, 0,
"3", "B", 1e1, 0,
"4", "A", 1e-3, 115,
"4", "B", 1e-3, 31,
"5", "A", 1e-3, 127,
"5", "B", 1e-4, 26,
"6", "A", 1e-3, 106,
"6", "B", 1e-2, 132,
"7", "A", 1e-3, 47,
"7", "B", 1e-2, 76,
"8", "A", 1e-3, 60,
"8", "B", 1e-3, 42)
cfuTable1B1 <- tribble(
~Number, ~BioRep, ~Dilution, ~Colonies,
"1", "A", 1e-4, 5,
"1", "B", 1e-5, 2,
"2", "A", 1e-1, 69,
"2", "B", 1e-1, 74,
"3", "A", 1e-4, 36,
"3", "B", 1e-4, 35,
"4", "A", 1e-3, 25,
"4", "B", 1e-4, 28,
"5", "A", 1e-3, 10,
"5", "B", 1e-4, 21,
"6", "A", 1e-3, 32,
"6", "B", 1e-3, 26,
"7", "A", 1e-3, 44,
"7", "B", 1e-3, 30,
"8", "A", 1e-4, 27,
"8", "B", 1e-4, 24)
cfuTable1B2 <- tribble(
~Number, ~BioRep, ~Dilution, ~Colonies,
"1", "A", 1e-3, 20,
"1", "B", 1e-3, 8,
"2", "A", 1e-3, 21,
"2", "B", 1e-3, 3,
"3", "A", 1e-4, 58,
"3", "B", 1e-4, 26,
"4", "A", 1e-4, 25,
"4", "B", 1e-2, 27,
"5", "A", 1e-3, 59,
"5", "B", 1e-3, 14,
"6", "A", 1e-3, 53,
"6", "B", 1e-3, 36,
"7", "A", 1e-3, 37,
"7", "B", 1e-3, 21,
"8", "A", 1e-1, 84,
"8", "B", 1e-3, 80,
"9", "A", 1e-2, 6,
"9", "B", 1e-1, 69)
cfuTable2A1 <- tribble(
~Number, ~BioRep, ~Dilution, ~Colonies,
"1", "A", 1e-4, 26,
"1", "B", 1e-4, 87,
"2", "A", 1e-4, 45,
"2", "B", 1e-5, 39,
"3", "A", NA, NA,
"3", "B", NA, NA)
cfuTable2A2 <- tribble(
~Number, ~BioRep, ~Dilution, ~Colonies,
"1", "A", 1e-5, 22,
"1", "B", 1e-5, 17,
"2", "A", 1e-4, 6,
"2", "B", 1e1, 0,
"3", "A", NA, NA,
"3", "B", NA, NA)
cfuTable3A1 <- tribble(
~Number, ~BioRep, ~Dilution, ~Colonies,
"1", "A", 1e-2, 66,
"1", "B", 1e-3, 12,
"2", "A", 1e-2, 73,
"2", "B", 1e-2, 95,
"3", "A", 1e-3, 33,
"3", "B", 1e-3, 22,
"4", "A", NA, NA,
"4", "B", NA, NA)
cfuTable3A2 <- tribble(
~Number, ~BioRep, ~Dilution, ~Colonies,
"1", "A", 1e-2, 66,
"1", "B", 1e-3, 12,
"2", "A", 1e-2, 73,
"2", "B", 1e-2, 95,
"3", "A", 1e-3, 33,
"3", "B", 1e-3, 22,
"4", "A", NA, NA,
"4", "B", NA, NA)
cfuTables <- list(cfuTable1A1, cfuTable1A2, cfuTable1B1, cfuTable1B2, cfuTable2A1, cfuTable2A2, cfuTable3A1, cfuTable3A2) %>%
map(~mutate(.x, CFU_mL=(Colonies)/(Dilution*0.007))) %>%
map2(., batches, ~mutate(.x, Batch=.y)) %>%
purrr::reduce(full_join, by = c("Number", "BioRep", "Dilution", "Colonies", "CFU_mL", "Batch"))
strainKeys <- midlog_raw_data %>% # get strain information
select(Number, Batch, Strain, Species) %>%
unique()
cfuTables %>%
left_join(strainKeys) %>%
mutate(Protocol=str_sub(Batch, 1, 1),
BatchNum=str_sub(Batch, 3,3),
logCFU=log10(CFU_mL),
Strain=factor(Strain, levels=strains),
Species=factor(Species, levels=species)) %>%
ggplot(aes(x=Strain, y=logCFU, color=Species, fill=Species, shape=BatchNum, linetype=BatchNum)) +
geom_point() +
stat_summary_bin(fun = "mean", geom = "crossbar", show.legend = FALSE) +
scale_y_continuous(breaks = seq(1,10,1), limits = c(1,10), labels=c("10", expression(10^2), expression(10^3), expression(10^4), expression(10^5), expression(10^6), expression(10^7), expression(10^8), expression(10^9), expression(10^10))) +
theme(axis.text.x = element_text(angle = 45, hjust=1),
panel.grid.minor.y = element_blank(),
legend.position = "bottom") +
guides(shape=FALSE) +
labs(x=NULL, y="CFU/mL")
sessionInfo()
## R version 4.4.0 (2024-04-24)
## Platform: x86_64-apple-darwin20
## Running under: macOS Ventura 13.6.6
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.4-x86_64/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.4-x86_64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.0
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## time zone: America/New_York
## tzcode source: internal
##
## attached base packages:
## [1] grid stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] gridGraphics_0.5-1 ggpubr_0.6.0 cowplot_1.1.3 ggbeeswarm_0.7.2
## [5] corrplot_0.92 formattable_0.2.1 ape_5.8 rstatix_0.7.2
## [9] broom_1.0.6 xlsx_0.6.5 lubridate_1.9.3 forcats_1.0.0
## [13] stringr_1.5.1 dplyr_1.1.4 purrr_1.0.2 readr_2.1.5
## [17] tidyr_1.3.1 tibble_3.2.1 ggplot2_3.5.1 tidyverse_2.0.0
##
## loaded via a namespace (and not attached):
## [1] gtable_0.3.5 beeswarm_0.4.0 xfun_0.46 bslib_0.7.0
## [5] htmlwidgets_1.6.4 rJava_1.0-11 lattice_0.22-6 tzdb_0.4.0
## [9] vctrs_0.6.5 tools_4.4.0 generics_0.1.3 parallel_4.4.0
## [13] fansi_1.0.6 highr_0.11 pkgconfig_2.0.3 lifecycle_1.0.4
## [17] farver_2.1.2 compiler_4.4.0 munsell_0.5.1 carData_3.0-5
## [21] vipor_0.4.7 htmltools_0.5.8.1 sass_0.4.9 yaml_2.3.9
## [25] crayon_1.5.3 pillar_1.9.0 car_3.1-2 jquerylib_0.1.4
## [29] cachem_1.1.0 abind_1.4-5 nlme_3.1-165 tidyselect_1.2.1
## [33] digest_0.6.36 stringi_1.8.4 labeling_0.4.3 fastmap_1.2.0
## [37] colorspace_2.1-0 cli_3.6.3 magrittr_2.0.3 xlsxjars_0.6.1
## [41] utf8_1.2.4 withr_3.0.0 scales_1.3.0 backports_1.5.0
## [45] bit64_4.0.5 timechange_0.3.0 rmarkdown_2.27 bit_4.0.5
## [49] ggsignif_0.6.4 zoo_1.8-12 hms_1.1.3 evaluate_0.24.0
## [53] knitr_1.48 rlang_1.1.4 Rcpp_1.0.13 glue_1.7.0
## [57] vroom_1.6.5 rstudioapi_0.16.0 jsonlite_1.8.8 R6_2.5.1